'''
 * @desc : 逻辑回归 梯度上升法
 * @auth : TYF
 * @date : 2019/9/13 - 15:58

  logistic回归:根据现有数据对分类边界线建立回归公式，以此进行分类

  如果x为分类器输入,w是权重,y是一条分类曲线
  y = wx + b

  将y使用sigmoid函数激活,得到一个0~1之间的数,但并不是我们想要的,我们想要二分类即应该得到二值化的结果,所以使用对数几率函数(以中sigmoid函数)
  y = 1 / ( 1 + exp(-(wx+b) ) )

  变形得到
  ln (y/1-y) = wx + b

  其中ln(y/1-y)就是对数几率,当y是正例概率那么1-y就是负例概率,则公式变形得到
  ln ( p(正) / p(负) ) = wx +b

  再变回去分别得到p(正)和p(负)
  p(y=1|x) = exp(wx+b) / ( 1 + exp(wx+b) )
  p(y=0|x) =     1     / ( 1 + exp(wx+b) )
  则合并成伯努利分布形式
  p(y|x) = p(y=1|x)^y + (1-p(y=1|x))^(1-y)

  接下来就是给定训练集求w,用均方差作为损失函数使用最小二乘法找一组w和b使得损失函数极小化,但是函数非凸有很多局部最小解,所以通常采用极大似然发求解
  极大似然法求w和b,已知结果反推导致结果最有可能的参数取值,取w为待估计参数θ,令样本x属于其真实类别的概率越大越好


'''

from __future__ import print_function
from numpy import *
import matplotlib.pyplot as plt


#TestSet.txt文本文件
#-0.017612	14.053064	0
#-1.395634	4.662541	1
#-0.752157	6.538620	0
#-1.322371	7.152853	0
#0.423363	11.054677	0

# 加载样本
def loadDataSet(file_name):
    # dataMat为原始数据,labelMat为原始数据的标签
    dataMat = []
    labelMat = []
    fr = open(file_name)
    for line in fr.readlines():
        lineArr = line.strip().split()
        if len(lineArr) == 1:
            continue    # 这里如果就一个空的元素，则跳过本次循环
        # 为了方便计算，我们将 X0 的值设为 1.0 ，也就是在每一行的开头添加一个 1.0 作为 X0
        dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])])
        labelMat.append(int(lineArr[2]))
    return dataMat, labelMat


# sigmoid跳跃函数
def sigmoid(inX):
    # return 1.0 / (1 + exp(-inX))
    # Tanh是Sigmoid的变形，与 sigmoid 不同的是，tanh 是0均值的。因此，实际应用中，tanh 会比 sigmoid 更好。
    return 2 * 1.0/(1+exp(-2*inX)) - 1


# 梯度上升
# dataMatIn 输入样本
# classLabels 样本标签
def gradAscent(dataMatIn, classLabels):
    # 样本转矩阵(一行一个样本)
    dataMatrix = mat(dataMatIn)  # 转换为 NumPy 矩阵
    # 标签转列向量
    labelMat = mat(classLabels).transpose()
    # m行(样本个数) n列(特征个数) 这里是5行3列
    m, n = shape(dataMatrix)
    # 梯度上升步长
    alpha = 0.0001
    # 迭代次数
    maxCycles = 5000
    # 初始化w 3行1列
    weights = ones((n, 1))
    # 开始迭代
    for k in range(maxCycles):
        # 计算预测值 5*3 乘 3*1 得到 5*1 并激活
        h = sigmoid(dataMatrix * weights)
        # 计算loss 5*1
        error = (labelMat - h)
        # 计算更新后的w
        # f(w)=wx对w求导等于x即输入样本
        # 所以调整量为 0.001 * 3*5 * 5*1是一个3*1向量
        # 这里发现每次更新w,样本矩阵都会参与计算效率低下
        weights = weights + alpha * dataMatrix.transpose() * error
    return array(weights)


# 随机梯度下降
def stocGradAscent0(dataMatrix, classLabels):
    # m行(样本个数) n列(特征个数) 这里是5行3列
    m, n = shape(dataMatrix)
    alpha = 0.01
    # 初始化w 三维数组
    weights = ones(n)
    # 开始迭代
    for i in range(m):
        # 随机梯度下降一次只用一个样本点来更新回归系数
        # 每次取第i个样本 1*3 * 3*1 得到一个具体数值
        h = sigmoid(sum(dataMatrix[i] * weights))
        # 计算该样本的loss 具体数值
        error = classLabels[i] - h
        # 0.01*(1*1)*(1*3) 得到三维数组
        weights = weights + alpha * error * dataMatrix[i]
    return weights


# 随机梯度下降
def stocGradAscent1(dataMatrix, classLabels, numIter=150):
    # 样本 5*3
    m, n = shape(dataMatrix)
    weights = ones(n)
    # 迭代150次观察是否收敛
    for j in range(numIter):
        # [0, 1, 2 .. m-1]
        dataIndex = range(m)
        for i in range(m):
            # 让alpha随着i和j的不断增大而减小(大于0.0001常数项)相当于学习到后期降低学习速率
            alpha = 4 / (
                1.0 + j + i
            ) + 0.0001
            # 随机产生一个 0～len()之间的一个值,即样本集中随机找一个样本
            randIndex = int(random.uniform(0, len(dataIndex)))
            # 随机样本 1*3 * 3*1 得到一个具体数值
            h = sigmoid(sum(dataMatrix[dataIndex[randIndex]] * weights))
            # 计算loss 得到一个具体数值
            error = classLabels[dataIndex[randIndex]] - h
            # 更新w
            weights = weights + alpha * error * dataMatrix[dataIndex[randIndex]]
            # 参与过更新w的样本从样本集中去掉
            # del (dataIndex[randIndex])
    return weights


# 可视化展示
def plotBestFit(dataArr, labelMat, weights):
    '''
        Desc:
            将我们得到的数据可视化展示出来
        Args:
            dataArr:样本数据的特征
            labelMat:样本数据的类别标签，即目标变量
            weights:回归系数
        Returns:
            None
    '''

    n = shape(dataArr)[0]
    xcord1 = []
    ycord1 = []
    xcord2 = []
    ycord2 = []
    for i in range(n):
        if int(labelMat[i]) == 1:
            xcord1.append(dataArr[i, 1])
            ycord1.append(dataArr[i, 2])
        else:
            xcord2.append(dataArr[i, 1])
            ycord2.append(dataArr[i, 2])
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.scatter(xcord1, ycord1, s=30, c='red', marker='s')
    ax.scatter(xcord2, ycord2, s=30, c='green')
    x = arange(-3.0, 3.0, 0.1)
    """
    y的由来，卧槽，是不是没看懂？
    首先理论上是这个样子的。
    dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])])
    w0*x0+w1*x1+w2*x2=f(x)
    x0最开始就设置为1叻， x2就是我们画图的y值，而f(x)被我们磨合误差给算到w0,w1,w2身上去了
    所以： w0+w1*x+w2*y=0 => y = (-w0-w1*x)/w2   
    """
    y = (-weights[0] - weights[1] * x) / weights[2]
    ax.plot(x, y)
    plt.xlabel('X')
    plt.ylabel('Y')
    plt.show()


def simpleTest():
    # 1.收集并准备数据
    dataMat, labelMat = loadDataSet("TestSet.txt")

    # print dataMat, '---\n', labelMat
    # 2.训练模型，  f(x)=a1*x1+b2*x2+..+nn*xn中 (a1,b2, .., nn).T的矩阵值
    # 因为数组没有是复制n份， array的乘法就是乘法
    dataArr = array(dataMat)
    # print dataArr
    # weights = gradAscent(dataArr, labelMat)
    # weights = stocGradAscent0(dataArr, labelMat)
    weights = stocGradAscent1(dataArr, labelMat)
    # print '*'*30, weights

    # 数据可视化
    plotBestFit(dataArr, labelMat, weights)



if __name__ == "__main__":
    simpleTest()
